using System;
using System.Windows.Forms;
using System.Drawing;
using System.Drawing.Drawing2D;
using System.IO;
using System.Runtime.Serialization;
using System.Runtime.Serialization.Formatters.Binary;
using System.Collections;
using System.Reflection;
using System.Security.Permissions;
using System.ComponentModel;
using System.Drawing.Imaging;
using System.Runtime.InteropServices;
using System.Data;
using System.Data.Odbc;
using NumericalRecipes;
using FftLib;
using BmgNet;
using StaticUtility;
using EngMATLib;
/*
namespace ObjectLib
{
	/// <summary>
	/// Uses a kernel analysis to evaluate the maximum accuracy possible for a given data set.
	/// </summary>
	[Serializable]
	public class KernAcc
	{
		/// <summary>
		/// Holds resources in memory.
		/// </summary>
		[NonSerialized]
		public static Nr NrStatic;
		/// <summary>
		/// A cache of lookup values for the kernel function.
		/// </summary>
		public static ExpSqrtLookup Lookup;
		/// <summary>
		/// The length of each ribbon.
		/// </summary>
		public int RibLen;
		/// <summary>
		/// An inventory of unique project numbers included in this analysis.
		/// </summary>
		public int[] ProjNums;
		/// <summary>
		/// The number of instances in each project (by row) and each ribbon type (by col, Boundary, Texture, or Shadow)
		/// </summary>
		[NonSerialized]
		public int[,] ProjCounts;
		/// <summary>
		/// The indexes for each project (by row) and each ribbon type (by col, Boundary, Texture, or Shadow)
		/// </summary>
		[NonSerialized]
		public int[,][] ProjIndices;
		/// <summary>
		/// A column for each type of information, a row for each ribbon that is a surface boundary.
		/// </summary>
		public double[,] BoundaryRibs;
		/// <summary>
		/// A column for each type of information, a row for each ribbon that is a surface texture.
		/// </summary>
		public double[,] TextureRibs;
		/// <summary>
		/// A column for each type of information, a row for each ribbon that is a shadow boundary.
		/// </summary>
		public double[,] ShadowRibs; 
		/// <summary>
		/// Matches each row of BoundaryRibs to an image (i.e. to an index of ProjNums)
		/// </summary>
		public int[] BoundaryProj;
		/// <summary>
		/// Matches each row of TextureRibs to an image (i.e. to an index of ProjNums)
		/// </summary>
		public int[] TextureProj;
		/// <summary>
		/// Matches each row of ShadowInfo to an image (i.e. to an index of ProjNums)
		/// </summary>
		public int[] ShadowProj;
		/// <summary>
		/// Static constructor holds numerical resources in memory.
		/// </summary>
		static KernAcc()
		{
			NrStatic = new Nr();
			//Lookup = new ExpSqrtLookup(0.01, -30.0, 0.0);
			Lookup = new ExpSqrtLookup(0.01, -30.0, 0.0);
		}
		/// <summary>
		/// Empty constructor
		/// </summary>
		public KernAcc() {}
		/// <summary>
		/// The number in each category (set in constructor);
		/// </summary>
		public int[] NperCat;
		/// <summary>
		/// Factors for correcting the likelihood calculation for base-rates so that priors are removed.
		/// </summary>
		public double[] BaseRateCorrection;
		/// <summary>
		/// Construct from information in a mat file.
		/// </summary>
		/// <param name="matFilename"></param>
		public static KernAcc FromMatFile(string matFilename)
		{
			KernAcc output = new KernAcc();
			int i;
			string fieldname;
			output.NperCat = new int[3];
			output.BaseRateCorrection = new double[3];
			int nTotal = 0;

			//	Open the mat file
			MATFile file = new MATFile(matFilename, EngMATLib.FileAccess.Read);

			//	Single number, the length of each ribbon.
			MatrixDescription mi = file.GetMatrixInfo("RibLen");
			double[,] temp = new double[1,1];
			file.GetMatrix("RibLen", ref temp);
			output.RibLen = (int)(temp[0,0]+0.5);

			//	Column vector of double where each double is a project number
			mi = file.GetMatrixInfo("ProjNums");
			temp = new double[mi.Rows, mi.Cols];
			file.GetMatrix("ProjNums", ref temp);
			output.ProjNums = new int[mi.Rows];
			for(i=0; i<mi.Rows; i++)
				output.ProjNums[i] = (int)(temp[i,0]+0.5);

			//	Columns of information where each row is a different ribbon
			fieldname = "Boundary";
			mi = file.GetMatrixInfo(fieldname);
			output.BoundaryRibs = new double[mi.Rows,mi.Cols];
			output.NperCat[0] = mi.Rows;
			file.GetMatrix(fieldname, ref output.BoundaryRibs);
			//	Matches each row of information to a particular image.
			fieldname = "BoundaryImg";
			mi = file.GetMatrixInfo(fieldname);
			temp = new double[mi.Rows,1];
			file.GetMatrix(fieldname, ref temp);
			output.BoundaryProj = new int[mi.Rows];
			for(i=0; i<mi.Rows; i++)
				output.BoundaryProj[i] = (int)(temp[i,0]-0.5);

			//	Columns of information where each row is a different ribbon
			fieldname = "Marking";
			mi = file.GetMatrixInfo(fieldname);
			output.TextureRibs = new double[mi.Rows,mi.Cols];
			output.NperCat[1] = mi.Rows;
			file.GetMatrix(fieldname, ref output.TextureRibs);
			//	Matches each row of information to a particular image.
			fieldname = "MarkingImg";
			mi = file.GetMatrixInfo(fieldname);
			temp = new double[mi.Rows,1];
			file.GetMatrix(fieldname, ref temp);
			output.TextureProj = new int[mi.Rows];
			for(i=0; i<mi.Rows; i++)
				output.TextureProj[i] = (int)(temp[i,0]-0.5);

			//	Columns of information where each row is a different ribbon
			fieldname = "Shadow";
			mi = file.GetMatrixInfo(fieldname);
			output.ShadowRibs = new double[mi.Rows,mi.Cols];
			output.NperCat[2] = mi.Rows;
			file.GetMatrix(fieldname, ref output.ShadowRibs);
			//	Matches each row of information to a particular image.
			fieldname = "ShadowImg";
			mi = file.GetMatrixInfo(fieldname);
			temp = new double[mi.Rows,1];
			file.GetMatrix(fieldname, ref temp);
			output.ShadowProj = new int[mi.Rows];
			for(i=0; i<mi.Rows; i++)
				output.ShadowProj[i] = (int)(temp[i,0]-0.5);

			nTotal = output.NperCat[0] + output.NperCat[1] + output.NperCat[2];
			double normer = 0.0;
			for(i=0; i<3; i++)
			{
				if( output.NperCat[i]==0 )
					output.BaseRateCorrection[i] = 0.0;
				else
				{
					output.BaseRateCorrection[i] = (double)nTotal/(double)output.NperCat[i];
					normer += output.BaseRateCorrection[i];
				}
			}
			output.BaseRateCorrection[0] /= normer;
			output.BaseRateCorrection[1] /= normer;
			output.BaseRateCorrection[2] /= normer;

			//	Return output
			return output;
		}
		/// <summary>
		/// Sets member variables "ProjCounts" and "ProjIndices".
		/// </summary>
		public void CountAndIndex()
		{
			int i,j;
			int[] ct = new int[ProjNums.Length];

			ProjCounts = new int[ProjNums.Length, 3];
			for(i=0; i<ProjNums.Length; i++)
				for(j=0; j<3; j++)
					ProjCounts[i,j] = 0;

			//	Count boundaries
			for(i=0; i<BoundaryProj.Length; i++)
				ProjCounts[BoundaryProj[i],0]++;
			//	Count textures
			for(i=0; i<TextureProj.Length; i++)
				ProjCounts[TextureProj[i],1]++;
			//	Count shadows
			for(i=0; i<ShadowProj.Length; i++)
				ProjCounts[ShadowProj[i],2]++;

			//	Allocate indices and evidence
			ProjIndices = new int[ProjNums.Length, 3][];
			for(j=0; j<ProjNums.Length; j++)
			{
				ProjIndices[j,0] = new int[ProjCounts[j,0]];
				ProjIndices[j,1] = new int[ProjCounts[j,1]];
				ProjIndices[j,2] = new int[ProjCounts[j,2]];
			}
			//	Index boundaries
			for(j=0; j<ProjNums.Length; j++)
				ct[j]=0;
			for(i=0; i<BoundaryProj.Length; i++)
				ProjIndices[BoundaryProj[i],0][ct[BoundaryProj[i]]++] = i;
			//	Index textures
			for(j=0; j<ProjNums.Length; j++)
				ct[j]=0;
			for(i=0; i<TextureProj.Length; i++)
				ProjIndices[TextureProj[i],1][ct[TextureProj[i]]++] = i;
			//	Index shadows
			for(j=0; j<ProjNums.Length; j++)
				ct[j]=0;
			for(i=0; i<ShadowProj.Length; i++)
				ProjIndices[ShadowProj[i],2][ct[ShadowProj[i]]++] = i;
		}
		/// <summary>
		/// Return a fit value for a particular project and particular standard deviation.
		/// </summary>
		/// <param name="projNum"></param>
		/// <param name="ksd"></param>
		/// <returns></returns>
		public double Fit(int projNum, double[] ksd)
		{
			return 0.0;
		}
		[NonSerialized]
		private double _dx, _dist, _wsum;
		[NonSerialized]
		private double[] _jlWorst;
		[NonSerialized]
		private int _pn, _cn, _ctProj, _ir, _dim, _co;
		[NonSerialized]
		private int[] _indices;
		[NonSerialized]
		private double[,] _catRibs=null;
		[NonSerialized]
		private double[] _weights, _zeros, _coordinate, _ksd;
		/// <summary>
		/// Computes the weights (towards each category) for a particular coordinate.  Samples from "projNum" are
		/// excluded from this calculation.  Assumes that CountAndIndex has already been called.
		/// </summary>
		/// <param name="projExclude">The project number to be excluded.  Set to -1 if all are to be included.</param>
		/// <param name="ksd">The kernel standard deviation for each dimension.</param>
		/// <param name="coordinate">A coordinate the same length as ksd.</param>
		/// <param name="weights">Accumulated evidence for membership in each category (Boundary, Texture, or Shadow).
		/// This should already be initialized to zeros.  This vector is not normalized.</param>
		public void Weights(int projExclude, double[] lnksd, double[] coordinate, double[] weights)
		{
			//	Get kernel variance
			for(_co=0; _co<lnksd.Length; _co++)
				_ksd[_co] = Math.Exp(lnksd[_co]);
			//	For each project
			for(_pn = 0; _pn < ProjNums.Length; _pn++)
			{
				//	That is not the project to be excluded
				if( _pn!= projExclude )
				{
					//	For each category
					for(_cn=0; _cn < 3; _cn++)
					{
						//	Get the relevant data
						if(_cn==0)
							_catRibs = this.BoundaryRibs;
						else if(_cn==1)
							_catRibs = this.TextureRibs;
						else
							_catRibs = this.ShadowRibs;

						//	For each ribbon index
						_ctProj = this.ProjCounts[_pn,_cn];
						_indices = ProjIndices[_pn,_cn];
						for(_ir=0; _ir < _ctProj; _ir++)
						{
							//	Compute distance squared
							_dist = 0.0;
							//	For each coordinate used
							for(_co=0; _co<lnksd.Length; _co++)
							{
								_dx = (coordinate[_co]-_catRibs[_indices[_ir],UsingColumns[_co]])/_ksd[_co];
								_dist+= _dx*_dx;
							}
							//	Compute Gaussian weight
							weights[_cn]+= Lookup.Eval(-_dist);
						}
					}
				}
			}
			if( this.CorrectingBaseRates )
			{
				weights[0] *= this.BaseRateCorrection[0];
				weights[1] *= this.BaseRateCorrection[1];
				weights[2] *= this.BaseRateCorrection[2];
			}
		}
		public double[] Ksd, InitialKsd;
		public void PrepareSolveKsd()
		{
			_dim = this.BoundaryRibs.GetLength(1);
			_weights = new double[3];
			_zeros = new double[_dim];
			_coordinate = new double[_dim];
			_ksd = new double[_dim];
			_jlWorst = new double[3];
			int i;

			InitialKsd = new double[_dim];
			double[] lnksd = new double[_dim];

			//	Boundary
			double[] temp = Util.Means(this.BoundaryRibs);
			temp = Util.Sds(this.BoundaryRibs, temp);
			for(i=0; i<_dim; i++)
				InitialKsd[i]+=temp[i];

			//	Texture
			temp = Util.Means(this.TextureRibs);
			temp = Util.Sds(this.TextureRibs, temp);
			for(i=0; i<_dim; i++)
				InitialKsd[i]+=temp[i];

			//	Shadow
			temp = Util.Means(this.ShadowRibs);
			temp = Util.Sds(this.ShadowRibs, temp);
			for(i=0; i<_dim; i++)
				InitialKsd[i]+=temp[i];

			for(i=0; i<_dim; i++)
				InitialKsd[i] /= 10.0;

			//Ksd = new double[] {0.110430106, 0.136501239, 0.137493144, 0.134196529, 0.057672173, 0.046531576, 0.030081232, 0.026448159, 0.07873665};
		}
		public int[] UsingColumns;
		public bool CorrectingBaseRates;
		public void SolveKsd(int[] UsingColumns, bool CorrectingBaseRates)
		{	
			_jlWorst[0] = _jlWorst[1] = _jlWorst[2] = 2.0*Math.Log(1.0/3.0);
			this.UsingColumns = UsingColumns;
			this.CorrectingBaseRates = CorrectingBaseRates;
			if( this.CorrectingBaseRates )
			{
				_jlWorst[0] = this.BaseRateCorrection[0]*Math.Log(1.0/3.0);
				_jlWorst[1] = this.BaseRateCorrection[1]*Math.Log(1.0/3.0);
				_jlWorst[2] = this.BaseRateCorrection[2]*Math.Log(1.0/3.0);
			}
			int i;
			double[] lnksd = new double[UsingColumns.Length];
			//	Normalize
			for(i=0; i<UsingColumns.Length; i++)
			{
				if( InitialKsd[UsingColumns[i]]>0.0 )
					lnksd[i] = Math.Log( InitialKsd[UsingColumns[i]] );
				else
					lnksd[i] = 0.0;
			}

			//	Allocate solver
			double chisqTol = (double)(this.NperCat[0]+this.NperCat[1]+this.NperCat[2])/1000.0;
			Im = new IngMin(chisqTol, lnksd, new IngMin.DelChisq(Fit));
			Im.Solve();

			//	Prepare output Ksd for each column in UsingColumns
			Ksd = new double[Im.Asol.Length];

			//	The kernel standard deviation is an Exp transform of the free parameters
			for(i=0; i<Im.Asol.Length; i++)
				Ksd[i] = Math.Exp(Im.Asol[i]);
		}
		[NonSerialized]
		public IngMin Im;
		/// <summary>
		/// Returns chi-square (-2*jointLikelihood) proportional to the amount of improperly accounted variance.
		/// </summary>
		/// <param name="ksd"></param>
		/// <returns></returns>
		public double Fit(double[] lnksd)
		{
			//	The joint likelihood
			double jl = 0.0;
			int pn, cn, ir, co, ctProj;
			double[] weights = new double[3];
			double[,] catRibs=null;
			int[] indices;
			//	For each project
			for(pn=0; pn < ProjNums.Length; pn++)
			{
				//	For each category
				for(cn=0; cn < 3; cn++)
				{
					//	Get the relevant data
					if(cn==0)
						catRibs = this.BoundaryRibs;
					else if(cn==1)
						catRibs = this.TextureRibs;
					else
						catRibs = this.ShadowRibs;

					//	For each ribbon index
					ctProj = this.ProjCounts[pn,cn];
					indices = this.ProjIndices[pn,cn];
					for(ir=0; ir < ctProj; ir++)
					{
						//	Read coordinates
						for(co=0; co<lnksd.Length; co++)
							_coordinate[co] = catRibs[indices[ir],UsingColumns[co]];
						//	Get weights
						weights[0] = weights[1] = weights[2] = 0.0;
						this.Weights(pn, lnksd, _coordinate, weights);
						// Normalize
						_wsum = weights[0] + weights[1] + weights[2];
						if(_wsum==0.0 )
						{
							jl+=_jlWorst[cn];
						}
						else
						{
							weights[0]/=_wsum;
							weights[1]/=_wsum;
							weights[2]/=_wsum;
							if( CorrectingBaseRates )
							{
								//	Sum the improperly accounted likelihood
								if( cn==0 )
								{
									jl+= this.BaseRateCorrection[0] *
										(weights[1]*Math.Log(Math.Max(1e-30,weights[1]))
										+weights[2]*Math.Log(Math.Max(1e-30,weights[2])));
								}
								else if( cn==1 )
								{
									jl+= this.BaseRateCorrection[1] *
										(weights[0]*Math.Log(Math.Max(1e-30,weights[0]))
										+weights[2]*Math.Log(Math.Max(1e-30,weights[2])));
								}
								else	//	cn==2
								{
									jl+= this.BaseRateCorrection[2] *
										(weights[1]*Math.Log(Math.Max(1e-30,weights[1]))
										+weights[0]*Math.Log(Math.Max(1e-30,weights[0])));
								}
							}
							else
							{
								//	Sum the improperly accounted likelihood
								if( cn==0 )
								{
									jl+= 
										weights[1]*Math.Log(Math.Max(1e-30,weights[1]))
										+weights[2]*Math.Log(Math.Max(1e-30,weights[2]));
								}
								else if( cn==1 )
								{
									jl+= 
										weights[0]*Math.Log(Math.Max(1e-30,weights[0]))
										+weights[2]*Math.Log(Math.Max(1e-30,weights[2]));
								}
								else	//	cn==2
								{
									jl+= 
										weights[1]*Math.Log(Math.Max(1e-30,weights[1]))
										+weights[0]*Math.Log(Math.Max(1e-30,weights[0]));
								}
							}
						}
					}
				}
			}
			//	Return chisq
			return -2.0*jl;
		}
		/// <summary>
		/// Return statistics on categorization performance.
		/// </summary>
		/// <param name="projNum">The project number.</param>
		/// <param name="ksd">The log of kernel standard deviation for each dimension.</param>
		/// <param name="confusion">A categorization confusion matrix.
		/// Each row represents the actual category and each column represents the predicted category.
		/// Entries along the diagonal are a count of the correct categorizations of boundary, texture, and shadows.  The other entries
		/// are miscategorization errors.</param>
		/// <param name="wmat">A weight matrix.  Like the confusion matrix, but this sums the weights instead of the discrete categorization assignment.</param>
		/// <param name="wmatFromNorm">A weight matrix where the weights from each exemplar are normalized to a sum of 1 before being
		/// added to wmatFromNorm.</param>
		/// <param name="nUncategorized">The number of stimuli from each category where categorization was impossible (kernel sd was too small).</param>
		/// <param name="jointLikelihood">Returns the joint likelihood in a confusion format.</param>
		public void RunStats(
			out int[][,] Confusion, out double[][,] Wmat, out double[][,] WmatFromNorm,
			out int[][] Nuncategorized, 
			out double[][,] JointLikelihood,
			out string textDescription)
		{
			int[,] tConfusion=new int[3,3];
			int[] tNuncategorized = new int[3];
			double[,] tWmat=new double[3,3], tWmatFromNorm=new double[3,3], tJointLikelihood=new double[3,3];
			int i, projNum;
			double[] lnksd = new double[this.UsingColumns.Length];
			//	Normalize
			for(i=0; i<lnksd.Length; i++)
			{
				if( Ksd[i]>0.0 )
					lnksd[i] = Math.Log(Ksd[i]);
				else
					lnksd[i] = 0.0;
			}

			textDescription = "";

			Confusion = new int[this.ProjNums.Length][,];
			Wmat = new double[this.ProjNums.Length][,];
			WmatFromNorm = new double[this.ProjNums.Length][,];
			JointLikelihood = new double[this.ProjNums.Length][,];
			Nuncategorized = new int[this.ProjNums.Length][];
			for(projNum = 0; projNum < this.ProjNums.Length; projNum++)
			{
				int[,] confusion = new int[3,3];
				double[,] wmat = new double[3,3];
				double[,] wmatFromNorm = new double[3,3];
				int[] nUncategorized = new int[3];
				double[,] jointLikelihood = new double[3,3];
				int cn, ir, co,j;
				double[] weights = new double[3];
				double wsum;
				double[,] catRibs=null;
				double[] coordinate = new double[_dim];
				int chosen=0, ctProj;

				_dim = this.BoundaryRibs.GetLength(1);
				_weights = new double[_dim];
				_zeros = new double[_dim];
				_coordinate = new double[_dim];
				_ksd = new double[_dim];

				//	For each category
				for(cn=0; cn < 3; cn++)
				{
					//	Get the relevant data
					if(cn==0)
						catRibs = this.BoundaryRibs;
					else if(cn==1)
						catRibs = this.TextureRibs;
					else
						catRibs = this.ShadowRibs;

					//	For each ribbon index
					ctProj = this.ProjCounts[projNum,cn];
					for(ir=0; ir < ctProj; ir++)
					{
						//	Read coordinates
						for(co=0; co<UsingColumns.Length; co++)
							coordinate[co] = catRibs[ir,UsingColumns[co]];
						//	Get weights
						weights[0] = weights[1] = weights[2] = 0.0;
						//	This function only takes 'UsingColumns' into account, which is what we want.
						this.Weights(projNum, lnksd, coordinate, weights);
						//	Increment the weight matrix
						wmat[cn,0] += weights[0];
						wmat[cn,1] += weights[1];
						wmat[cn,2] += weights[2];
						//	Normalize the weight matrix
						wsum = weights[0] + weights[1] + weights[2];
						if( wsum>0.0 )
						{
							weights[0]/=wsum;
							weights[1]/=wsum;
							weights[2]/=wsum;
							//	What is the category chosen?
							if( weights[0] > weights[1] )
							{
								if( weights[0] > weights[2] )	//	category 0
									chosen = 0;
								else // category 2 
									chosen = 2;
							}
							else
							{
								if( weights[1] > weights[2] )	//	category 1
									chosen = 1;
								else //category 2
									chosen = 2;
							}
							confusion[cn,chosen]++;
							//	Increment the weight matrix from normalized weights
							wmatFromNorm[cn,0] += weights[0];
							wmatFromNorm[cn,1] += weights[1];
							wmatFromNorm[cn,2] += weights[2];
							//	Increment the joint likelihood
							jointLikelihood[cn,0] += weights[0]*Math.Log(Math.Max(1e-30,weights[0]));
							jointLikelihood[cn,1] += weights[1]*Math.Log(Math.Max(1e-30,weights[1]));
							jointLikelihood[cn,2] += weights[2]*Math.Log(Math.Max(1e-30,weights[2]));
						}
						else
						{
							weights[0] = weights[1] = weights[2] = 0.0;
							nUncategorized[cn]++;
							jointLikelihood[cn,0] += _jlWorst[cn];
							jointLikelihood[cn,1] += _jlWorst[cn];
							jointLikelihood[cn,2] += _jlWorst[cn];
						}
					}
				}
				for(i=0; i<3; i++)
				{
					tNuncategorized[i] += nUncategorized[i];
					for(j=0; j<3; j++)
					{
						tConfusion[i,j] += confusion[i,j];
						tWmat[i,j] += wmat[i,j];
						tWmatFromNorm[i,j] += wmatFromNorm[i,j];
						tJointLikelihood[i,j] += jointLikelihood[i,j];
					}
				}
				textDescription = textDescription + "\n------------------------------------------\n(" + projNum.ToString() + ")\nN-Uncategorized:\n";
				for(i=0; i<nUncategorized.Length; i++)
					textDescription = textDescription + nUncategorized[i].ToString() + "  ";
				textDescription = textDescription + "\nConfusion Matrix:\n" + Util.MatrixToString(confusion);
				textDescription = textDescription + "\nTotal Weight Matrix:\n" + Util.MatrixToString(wmat);
				textDescription = textDescription + "\nTotal Normalized Matrix:\n" + Util.MatrixToString(wmatFromNorm);
				textDescription = textDescription + "\nTotal Joint Likelihood:\n" + Util.MatrixToString(jointLikelihood) +" \n";
				Confusion[i] = confusion;
				Wmat[i] = wmat;
				WmatFromNorm[i] = wmatFromNorm;
				JointLikelihood[i] = jointLikelihood;
				Nuncategorized[i] = nUncategorized;
			}

			string prim = "\nSUMMARY============================================================================\nKernel Sd:\n";
			for(i=0; i<this.Ksd.Length; i++)
				prim = prim + Math.Exp(lnksd[i]).ToString() + "\n";
			prim = prim + '\n' + UsingColumns.Length.ToString() + " Columns Used, Column Numbers: ";
			for(i=0; i<this.UsingColumns.Length; i++)
				if(i== this.UsingColumns.Length-1)
					prim = prim + UsingColumns[i].ToString() + ".\n";
				else
					prim = prim + UsingColumns[i].ToString() + ", ";
			if( this.CorrectingBaseRates )
				prim = prim + "Corrected For Base-Rates (priors removed)";
			else
				prim = prim + "NOT Corrected For Base-Rates (priors included)";
			for(i=0; i<tNuncategorized.Length; i++)
				prim = prim + tNuncategorized[i].ToString() + "  ";
			prim = prim + "\nConfusion Matrix:\n" + Util.MatrixToString(tConfusion);
			prim = prim + "\nTotal Weight Matrix:\n" + Util.MatrixToString(tWmat);
			prim = prim + "\nTotal Normalized Matrix:\n" + Util.MatrixToString(tWmatFromNorm);
			prim = prim + "\nTotal Joint Likelihood:\n" + Util.MatrixToString(tJointLikelihood) +" \n";
			textDescription = prim + "\n\n+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n\n" + textDescription;
		}
	}
	public class BracketSorter : IComparer
	{
		public BracketSorter()	{}
		public int[] xx, yy;
		public int Compare(object x, object y)
		{
			xx = (int[])x;
			yy = (int[])y;
			if(xx[3]<yy[3])
				return -1;
			else if(xx[3]==yy[3])
				return 0;
			return 1;
		}
	}
	[Serializable]
	public class Sample : IComparable
	{
		public static int CompareDimension=0;
		public static Sample _compareTo;
		public static IComparer Sorter = new BracketSorter();
		/// <summary>
		/// The absolute maximum normalized distance for bracketing.
		/// </summary>
		public const double DIST_MIN = 5.0;
		public const double SQ_DIST_MIN = DIST_MIN * DIST_MIN;
		/// <summary>
		/// A cache of lookup values for the kernel function.
		/// </summary>
		public static ExpLookup Lookup = new ExpLookup(0.01, -SQ_DIST_MIN, 0.0);
		public int Tag;
		public double[] Location;
		public int[] ListIndices;
		public int Dim 
		{
			get
			{
				return Location.Length;
			}
		}
		public Sample(double[] Location, int Tag)
		{
			this.Location = Location;
			this.ListIndices = new int[Location.Length];
			this.Tag = Tag;
		}
		public int CompareTo(object obj)
		{
			_compareTo = (Sample)obj;
			if(Location[CompareDimension] < _compareTo.Location[CompareDimension])
				return -1;
			else if(Location[CompareDimension] > _compareTo.Location[CompareDimension])
				return 1;
			return 0;
		}
		/// <summary>
		/// </summary>
		/// <param name="distThresh">The distance threshold for each dimension.</param>
		/// <param name="output">Outputs the indices of the brackets that sit on the bounds of the range.
		/// output[.][0] is the lower bound of the bracket, output[.][1] is the upper bound, 
		/// output[.][2] is the bracket dimension (for compatibility with BracketSorter), and
		/// output[.][3] is the bracket range (for compatiblity with BracketSorter).</param>
		/// <param name="sort">If true, the brackets are sorted by a BracketSorter object.</param>
		public void FindBrackets(Distribution category, double[] distThresh, int[][] output, bool sort)
		{
			//	Surround the bracket
			double limUp, limDown, limTry, a, b, lim;
			int iUp, iDown, iTry, i;
			bool searchToggle = true;
			int[] br;
			Sample[] list;
			for(i=0; i<distThresh.Length; i++)
			{
				list = category.SortedSamples[i];
				br = output[i];
				lim = Location[i]-distThresh[i];

				//	Bracket downward
				iDown = 0;
				limDown = list[iDown].Location[i];
				if( limDown < lim )
				{
					iUp = list.Length-1;
					limUp = list[iUp].Location[i];
					if( lim<limDown )
						lim = limDown;
					else if(lim>limUp)
						lim = limUp;
					while( iDown<(iUp-1) && limDown < lim )
					{
						if( searchToggle )
						{
							a = lim-limDown;
							b = limUp-lim;
							iTry = iDown + (int)(a/(a+b)*(double)(iUp-iDown)+0.5);
						}
						else
							iTry = iDown + (int)(0.5*(double)(iUp-iDown)+0.5);
						searchToggle = !searchToggle;
						if( iTry==iDown )
							iTry = iDown + 1;
						else if( iTry==iUp )
							iTry = iUp - 1;
						limTry = list[iTry].Location[i];
						if(limTry >= lim)
						{
							iUp = iTry;
							limUp = limTry;
						}
						else if(limTry < lim)
						{
							iDown = iTry;
							limDown = limTry;
						}
					}
					br[0] = iUp;
				}
				else
					br[0] = iDown;
				
				//	Bracket upward
				lim = Location[i]+distThresh[i];
				iUp = list.Length-1;
				limUp = list[iUp].Location[i];
				if( limUp > lim )
				{
					iDown = br[0];
					if( iDown<list.Length )
					{
						limDown = list[iDown].Location[i];
						if( lim<limDown )
							lim = limDown;
						else if(lim>limUp)
							lim = limUp;
						while( iDown<(iUp-1) && limUp > lim )
						{
							if( searchToggle )
							{
								a = lim-limDown;
								b = limUp-lim;
								iTry = iDown + (int)(a/(a+b)*(double)(iUp-iDown)+0.5);
							}
							else
								iTry = iDown + (int)(0.5*(double)(iUp-iDown)+0.5);
							searchToggle = !searchToggle;
							if( iTry==iDown )
								iTry = iDown + 1;
							else if( iTry==iUp )
								iTry = iUp - 1;
							limTry = list[iTry].Location[i];
							if(limTry >= lim)
							{
								iUp = iTry;
								limUp = limTry;
							}
							else if(limTry < lim)
							{
								iDown = iTry;
								limDown = limTry;
							}
						}
						br[1] = iDown;
					}
					else
						br[1] = iUp;
				}
				else
					br[1] = iUp;

				//	Update the range and index for compatibility with BracketSorter
				br[2] = i;
				br[3] = br[1]-br[0]+1;
			}
			if(sort)
				Array.Sort( output, 0, distThresh.Length, Sorter );
		}
		/// <summary>
		/// The objects weight will be computed based on all samples inside the sorted brackets that are less than
		/// MIN_DIST away.
		/// </summary>
		/// <param name="category"></param>
		/// <param name="sd"></param>
		/// <param name="sortedBrackets"></param>
		/// <returns></returns>
		public double Weight(Distribution category, double[] sd, int[][] sortedBrackets)
		{
			Sample sam;
			Sample[] sams;
			double sqdist, dx, weight = 0.0;
			int i,j,dim;

			if( sortedBrackets[0][3] <= 0 )
				return 0.0;

			//	Use the first bracket to index the possibilities.
			bool withinRange;
			int[] br, brFirst = sortedBrackets[0];
			sams = category.SortedSamples[brFirst[2]];
			for(i=brFirst[0]; i<=brFirst[1]; i++)
			{
				sam = sams[i];
				dx = (sam.Location[brFirst[2]]-this.Location[brFirst[2]])/sd[brFirst[2]];
				sqdist = dx*dx;
				withinRange = true;
				for(dim = 1; dim<sd.Length; dim++)
				{
					//	The next bracket
					br = sortedBrackets[dim];
					//	The next dimension
					j = br[2];
					if( sam.ListIndices[j]>=br[0] && sam.ListIndices[j]<=br[1] )
					{
						dx = (sam.Location[j]-this.Location[j])/sd[j];
						sqdist += dx*dx;
						if( sqdist > SQ_DIST_MIN )
						{
							withinRange = false;
							break;
						}
					}
					else
					{
						withinRange = false;
						break;
					}
				}
				if( withinRange )
					weight += Lookup.Eval(-sqdist);
			}
			return weight;
		}
	}
	[Serializable]
	public class Distribution : ArrayList
	{
		public int Dim;
		public Distribution(int Dim) : base()
		{
			this.Dim = Dim;
		}
		public Distribution(int Dim, int Count) : base(Count)
		{
			this.Dim = Dim;
		}
		public Distribution(double[,] samples, int[] tags) : base(samples.GetLength(0))
		{
			int i,j;
			int len = samples.GetLength(0);
			this.Dim = samples.GetLength(1);
			double[] loc;
			for(i=0; i<len; i++)
			{
				loc = new double[Dim];
				for(j=0; j<Dim; j++)
					loc[j] = samples[i,j];
				this.Add( new Sample( loc, tags[i] ) );
			}
		}
		public new Sample this[int ind]
		{
			get
			{
				return (Sample)base[ind];
			}
			set
			{
				base[ind] = value;
			}
		}
		public Sample[][] SortedSamples;
		public void CalcSortedSamples()
		{
			int i,j;
			SortedSamples = new Sample[Dim][];
			Sample[] sam;
			for(i=0; i<Dim; i++)
			{
				Sample.CompareDimension = i;
				base.Sort();
				SortedSamples[i] = sam = (Sample[])this.ToArray(typeof(Sample));
				for(j=0; j<sam.Length; j++)
					sam[j].ListIndices[i] = j;
			}
		}
		/// <summary>
		/// Creates a copy of this distribution but with samples 'Tag' removed.  New Sample objects are created,
		/// but the Sample member variable 'Location' points to the same memory location as the original Sample objects.
		/// </summary>
		/// <param name="val"></param>
		/// <returns></returns>
		public void TagSeparated(int tag, out Distribution untagged, out Distribution tagged)
		{
			int t;
			untagged = new Distribution(Dim);
			tagged = new Distribution(Dim);
			for(int i=0; i<this.Count; i++)
			{
				t = this[i].Tag;
				if( t != tag )
					untagged.Add(new Sample(this[i].Location, t));
				else
					tagged.Add(new Sample(this[i].Location, t));
			}
		}
	}
	/// <summary>
	/// Holds a set of distributions separated by category and sample tags.
	/// Each untagged set provides a separate sorting of its samples for each dimension.
	/// This object facilitates the calculation of categorization accuracy for a multivariate, multi-category distribution;
	/// </summary>
	[Serializable]
	public class DistributionSet
	{
		public double[] Var, Ksd;
		public int Dim;
		public int Ncats;
		public int Ntags;
		public int[] NperCat;
		public int Ntotal;
		public double[] BaseRateCorrection;
		/// <summary>
		/// Distribution objects by [Category, Tag]
		/// </summary>
		public Distribution[,] Tagged;
		/// <summary>
		/// Distribution objects by [Category, Tag]
		/// </summary>
		public Distribution[,] Untagged;
		public DistributionSet( Distribution[] primary, int[] tags )
		{
			int i,j,d;
			Distribution p;
			Sample sam;
			Dim = primary[0].Dim;
			Var = new double[Dim];
			Ksd = new double[Dim];
			double[] mean;
			double dx;
			Ncats = primary.Length;
			Ntags = tags.Length;
			NperCat = new int[Ncats];
			BaseRateCorrection = new double[Ncats];
			Ntotal = 0;
			for(i=0; i<Ncats; i++)
			{
				p = primary[i];
				NperCat[i] = p.Count;
				Ntotal += NperCat[i];
				mean = new double[Dim];
				for(j=0; j<p.Count; j++)
				{
					sam = p[j];
					for(d=0; d<Dim; d++)
						mean[d] += sam.Location[d];
				}
				for(d=0; d<Dim; d++)
					mean[d] /= (double)NperCat[i];
				for(j=0; j<p.Count; j++)
				{
					sam = p[j];
					for(d=0; d<Dim; d++)
					{
						dx = sam.Location[d]-mean[d];
						Var[d] += dx*dx;
					}
				}
			}
			for(d=0; d<Dim; d++)
				Var[d] /= (double)(Ntotal-1);
			this.Tagged = new Distribution[Ncats,Ntags];
			this.Untagged = new Distribution[Ncats,Ntags];
			for(i=0; i<Ncats; i++)
			{
				for(j=0; j<Ntags; j++)
				{
					primary[i].TagSeparated(tags[j], out Untagged[i,j], out Tagged[i,j]);
					Untagged[i,j].CalcSortedSamples();
				}
			}
			
			//	Compute the base-rate correction factor (for removing priors)
			double normer = 0.0;
			for(i=0; i<Ncats; i++)
			{
				if( NperCat[i]==0 )
					BaseRateCorrection[i] = 0.0;
				else
				{
					BaseRateCorrection[i] = (double)Ntotal/(double)NperCat[i];
					normer += BaseRateCorrection[i];
				}
			}
			for(i=0; i<Ncats; i++)
				BaseRateCorrection[i] /= normer;
			distThresh = new double[Dim];
			brackets = new int[Dim][];
			for(i=0; i<Dim; i++)
				brackets[i] = new int[4];
		}
		[NonSerialized]
		public IngMin Im;
		public bool RemovePriors;
		private double[] _jlWorst, distThresh;
		private int[][] brackets;
		/// <summary>
		/// Returns chi-square (-2*jointLikelihood) proportional to the amount of improperly accounted variance.
		/// </summary>
		/// <param name="ksd"></param>
		/// <returns></returns>
		public double Fit(double[] lnksd)
		{
			int i,t,c,cc,nt;
			double p;
			Distribution dist, undist;
			Sample sam;
			distThresh = new double[lnksd.Length];
			for(i=0; i<lnksd.Length; i++)
			{
				Ksd[i] = Math.Exp(lnksd[i]);
				distThresh[i] = Sample.DIST_MIN*Ksd[i];
			}
			//	The joint likelihood
			double jl = 0.0;
			double[] weights;
			double wsum;
			//	For each category and tag
			for(c=0; c<Ncats; c++)
			{
				for(t=0; t<Ntags; t++)
				{
					dist = this.Tagged[c,t];
					nt = dist.Count;
					for(i=0; i<nt; i++)
					{
						sam = dist[i];
						weights = new double[Ncats];
						wsum = 0.0;
						for(cc=0; cc<Ncats; cc++)
						{
							undist = this.Untagged[cc,t];
							sam.FindBrackets(undist, distThresh, brackets, true);
							if( this.RemovePriors )
								weights[cc] = BaseRateCorrection[cc] * sam.Weight(undist, Ksd, brackets);
							else
								weights[cc] = sam.Weight(undist, Ksd, brackets);
							wsum += weights[cc];
						}
						if( wsum > 0.0 )
						{
							if( this.RemovePriors )
							{
								p = Math.Min(0.9999999999999999,Math.Max(0.0000000000000001,  weights[c]/wsum));
								jl += BaseRateCorrection[c] * Math.Log(p);
							}
							else
							{
								p = Math.Min(0.9999999999999999,Math.Max(0.0000000000000001,  weights[c]/wsum));
								jl += Math.Log(p);
							}
						}
						else
							jl += _jlWorst[c];
					}
				}
			}
			//	Return chisq
			return -2.0*jl;
		}
		public void SolveKsd(bool RemovePriors)
		{	
			this.RemovePriors = RemovePriors;

			_jlWorst = new double[Ncats];
			_jlWorst[0] = Math.Log(1.0/(double)Ncats);
			int i;
			for(i=1; i<Ncats; i++)
				_jlWorst[i] = _jlWorst[0];

			if( RemovePriors )
				for(i=0; i<Ncats; i++)
					_jlWorst[i]*=this.BaseRateCorrection[i];
			double[] lnksd = new double[Dim];
			//	Normalize
			for(i=0; i<Dim; i++)
			{
				if( Var[i]>0.0 )
				{
					Ksd[i] = Math.Sqrt(Var[i])/10.0 * (double)Dim;
					lnksd[i] = Math.Log( Ksd[i] );
				}
				else
				{
					Ksd[i] = 1.0;
					lnksd[i] = 0.0;
				}
			}

			//	Allocate solver
			double chisqTol = (double)(this.Ntotal)/500.0;
			Im = new IngMin(chisqTol, lnksd, new IngMin.DelChisq(Fit));
			Im.Solve();

			//	The kernel standard deviation is an Exp transform of the free parameters
			for(i=0; i<Im.Asol.Length; i++)
				Ksd[i] = Math.Exp(Im.Asol[i]);
		}
		/// <summary>
		/// Text statistics on categorization performance.
		/// </summary>
		public string Stats
		{
			get
			{
				int i,t,c,cc,nt,cmax;
				double wmax=0.0;
				Distribution dist, undist;
				Sample sam;
				distThresh = new double[Ksd.Length];
				for(i=0; i<Ksd.Length; i++)
					distThresh[i] = Sample.DIST_MIN*Ksd[i];

				int[,] tConfusion=new int[Ncats,Ncats];
				int[] tNuncategorized = new int[Ncats];
				for(c=0; c<Ncats; c++)
				{
					tNuncategorized[c] = 0;
					for(cc=0; cc<Ncats; cc++)
						tConfusion[c,cc] = 0;
				}

				//	The joint likelihood
				double[] weights;
				double wsum;
				//	For each category and tag
				for(c=0; c<Ncats; c++)
				{
					for(t=0; t<Ntags; t++)
					{
						dist = this.Tagged[c,t];
						nt = dist.Count;
						for(i=0; i<nt; i++)
						{
							sam = dist[i];
							weights = new double[Ncats];
							wsum = 0.0;
							cmax = 0;
							for(cc=0; cc<Ncats; cc++)
							{
								undist = this.Untagged[cc,t];
								sam.FindBrackets(undist, distThresh, brackets, true);
								if( this.RemovePriors )
									weights[cc] = BaseRateCorrection[cc] * sam.Weight(undist, Ksd, brackets);
								else
									weights[cc] = sam.Weight(undist, Ksd, brackets);
								wsum += weights[cc];
								if( cc==0 )
									wmax = weights[0];
								else
								{
									if(weights[cc] > wmax)
									{
										cmax = cc;
										wmax = weights[cc];
									}
								}
							}
							if(wmax==0.0)
								tNuncategorized[c]++;
							else
								tConfusion[c,cmax]++;
						}
					}
				}
				
				string txt = "Ksd:\n";
				for(i=0; i<Ksd.Length; i++)
					txt = txt + Ksd[i].ToString() + '\n';
				txt = txt + "\nConfusion:\n" + Util.MatrixToString(tConfusion);
				txt = txt + "\nUncategorized:\n";
				for(i=0; i<Ncats; i++)
					txt = txt + tNuncategorized[i].ToString() + '\n';

				if(this.Im!=null)
				{
					txt = txt + "\nSolver Output:\n" + Im.FullVerbose;
				}
				return txt;
			}
		}
		/// <summary>
		/// Construct from information in a mat file.
		/// </summary>
		/// <param name="matFilename"></param>
		public static DistributionSet FromMatFile(string matFilename, int maxParticles, int[] columns)
		{
			Distribution[] primary = new Distribution[3];
			double[,] ribs=null;
			int[] ribTags=null;
			int[] rows = null;
			double[,] temp=null;
			float ftemp;
			int i;
			string fieldname;

			//	Open the mat file
			MATFile file = new MATFile(matFilename, EngMATLib.FileAccess.Read);

			//	Column vector of associations to unique projects numbered 1 .. length(ProjNums)
			EngMATLib.MatrixDescription mi = file.GetMatrixInfo("ProjNums");
			int[] tags = new int[mi.Rows];
			for(i=0; i<mi.Rows; i++)
				tags[i] = (i+1);

			//	Columns of information where each row is a different ribbon
			fieldname = "Boundary";
			mi = file.GetMatrixInfo(fieldname);
			if( mi.Rows <= maxParticles )
			{
				rows = new int[mi.Rows];
				for(i=0; i<mi.Rows; i++)
					rows[i] = i;
			}
			else
			{
				rows = new int[maxParticles];
				ftemp = (float)mi.Rows / (float)maxParticles;
				for(i=0; i<maxParticles; i++)
					rows[i] = (int)((float)i*ftemp);
			}
			ribs = new double[rows.Length,columns.Length];
			ribTags = new int[rows.Length];
			file.GetMatrix(fieldname, ref ribs, rows, columns);
			//	Matches each row of information to a particular image.
			fieldname = "BoundaryImg";
			temp = new double[mi.Rows,1];
			file.GetMatrix(fieldname, ref temp, rows, new int[] {0});
			for(i=0; i<rows.Length; i++)
				ribTags[i] = (int)(temp[i,0]+0.5);			
			primary[0] = new Distribution(ribs, ribTags);

			//	Columns of information where each row is a different ribbon
			fieldname = "Marking";
			mi = file.GetMatrixInfo(fieldname);
			if( mi.Rows <= maxParticles )
			{
				rows = new int[mi.Rows];
				for(i=0; i<mi.Rows; i++)
					rows[i] = i;
			}
			else
			{
				rows = new int[maxParticles];
				ftemp = (float)mi.Rows / (float)maxParticles;
				for(i=0; i<maxParticles; i++)
					rows[i] = (int)((float)i*ftemp);
			}
			ribs = new double[rows.Length,columns.Length];
			ribTags = new int[rows.Length];
			file.GetMatrix(fieldname, ref ribs, rows, columns);
			//	Matches each row of information to a particular image.
			fieldname = fieldname + "Img";
			temp = new double[mi.Rows,1];
			file.GetMatrix(fieldname, ref temp, rows, new int[] {0});
			for(i=0; i<rows.Length; i++)
				ribTags[i] = (int)(temp[i,0]+0.5);			
			primary[1] = new Distribution(ribs, ribTags);

			//	Columns of information where each row is a different ribbon
			fieldname = "Shadow";
			mi = file.GetMatrixInfo(fieldname);
			if( mi.Rows <= maxParticles )
			{
				rows = new int[mi.Rows];
				for(i=0; i<mi.Rows; i++)
					rows[i] = i;
			}
			else
			{
				rows = new int[maxParticles];
				ftemp = (float)mi.Rows / (float)maxParticles;
				for(i=0; i<maxParticles; i++)
					rows[i] = (int)((float)i*ftemp);
			}
			ribs = new double[rows.Length,columns.Length];
			ribTags = new int[rows.Length];
			file.GetMatrix(fieldname, ref ribs, rows, columns);
			//	Matches each row of information to a particular image.
			fieldname = fieldname + "Img";
			temp = new double[mi.Rows,1];
			file.GetMatrix(fieldname, ref temp, rows, new int[] {0});
			for(i=0; i<rows.Length; i++)
				ribTags[i] = (int)(temp[i,0]+0.5);			
			primary[2] = new Distribution(ribs, ribTags);

			//	Create output and return.
			return new DistributionSet(primary, tags);
		}
		/// <summary>
		/// Construct from information in a mat file.
		/// </summary>
		/// <param name="matFilename"></param>
		public static DistributionSet FromMatFile2(string matFilename, int maxParticles, int[] columns)
		{
			Distribution[] primary = new Distribution[2];
			double[,] ribs=null;
			int[] ribTags=null;
			int[] rows = null;
			double[,] temp=null;
			float ftemp;
			int i;
			string fieldname;

			//	Open the mat file
			MATFile file = new MATFile(matFilename, EngMATLib.FileAccess.Read);

			//	Column vector of associations to unique projects numbered 1 .. length(ProjNums)
			EngMATLib.MatrixDescription mi = file.GetMatrixInfo("ProjNums");
			int[] tags = new int[mi.Rows];
			for(i=0; i<mi.Rows; i++)
				tags[i] = (i+1);

			//	Columns of information where each row is a different ribbon
			fieldname = "Boundary";
			mi = file.GetMatrixInfo(fieldname);
			if( mi.Rows <= maxParticles )
			{
				rows = new int[mi.Rows];
				for(i=0; i<mi.Rows; i++)
					rows[i] = i;
			}
			else
			{
				rows = new int[maxParticles];
				ftemp = (float)mi.Rows / (float)maxParticles;
				for(i=0; i<maxParticles; i++)
					rows[i] = (int)((float)i*ftemp);
			}
			ribs = new double[rows.Length,columns.Length];
			ribTags = new int[rows.Length];
			file.GetMatrix(fieldname, ref ribs, rows, columns);
			//	Matches each row of information to a particular image.
			fieldname = fieldname + "Img";
			temp = new double[mi.Rows,1];
			file.GetMatrix(fieldname, ref temp, rows, new int[] {0});
			for(i=0; i<rows.Length; i++)
				ribTags[i] = (int)(temp[i,0]+0.5);			
			primary[0] = new Distribution(ribs, ribTags);

			//	Columns of information where each row is a different ribbon
			fieldname = "Contour";
			mi = file.GetMatrixInfo(fieldname);
			if( mi.Rows <= maxParticles )
			{
				rows = new int[mi.Rows];
				for(i=0; i<mi.Rows; i++)
					rows[i] = i;
			}
			else
			{
				rows = new int[maxParticles];
				ftemp = (float)mi.Rows / (float)maxParticles;
				for(i=0; i<maxParticles; i++)
					rows[i] = (int)((float)i*ftemp);
			}
			ribs = new double[rows.Length,columns.Length];
			ribTags = new int[rows.Length];
			file.GetMatrix(fieldname, ref ribs, rows, columns);
			//	Matches each row of information to a particular image.
			fieldname = fieldname + "Img";
			temp = new double[mi.Rows,1];
			file.GetMatrix(fieldname, ref temp, rows, new int[] {0});
			for(i=0; i<rows.Length; i++)
				ribTags[i] = (int)(temp[i,0]+0.5);			
			primary[1] = new Distribution(ribs, ribTags);

			//	Create output and return.
			return new DistributionSet(primary, tags);
		}
	}
}
*/